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The  fracture  of  Si  particles  due  to  internal  stresses  formed  during  the  intercalation  of  lithium  ions  was 
described  by  means  of  a  thermal  analogy  model  and  brittle  fracture  damage  parameter.  The  stresses  were 
calculated  following  the  diffusion  equation  and  equations  of  elasticity  with  an  appropriate  volumetric 
expansion  term.  The  results  were  compared  with  the  acoustic  emission  data  from  the  experiments  on 
electrochemical  cycling  of  Li  ion  half-cells  with  silicon  electrodes.  A  good  correlation  between  experiment 
and  prediction  was  observed.  The  results  of  computations  with  different  particle  sizes  show  the  existence 
of  a  critical  size  below  which  fracture  during  the  lithiation  is  not  expected  to  occur.  Such  a  critical  size 
appears  to  be  within  micrometer  scale. 
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1.  Introduction 

Electrification  of  automotive  drivetrains  or  powertrains  is  one 
of  the  primary  research  directions  for  the  reduction  of  petroleum 
dependence.  Lithium  ion  based  electrochemical  energy  storage 
technology  is  the  primary  candidate  for  plug-in  hybrid  electric 
vehicles  (PHEV)  and  fully  electric  vehicles  (EV),  and  an  ongoing  col¬ 
laborative  research  between  companies  and  research  institutions 
provides  evidence  for  a  high  potential  of  application  of  Li-ion  bat¬ 
teries  in  vehicle  development.  Despite  attractive  features  and  wide 
range  of  applications  in  smaller  scale  devices  such  as  consumer 
electronics,  the  service  life  of  Li-ion  battery  systems  is  still  lim¬ 
ited  due  to  the  processes  of  degradation  of  electrode  material  with 
repeated  charging  and  discharging,  and  mechanisms  of  such  degra¬ 
dation  have  not  been  sufficiently  investigated  and  understood. 

The  working  principle  of  Li  ion  battery  is  based  on 
repeated  transport  of  lithium  ions  through  electrolyte  and  inser¬ 
tion/extraction  from  the  electrode  host  structure  (intercalation 
compound).  It  has  been  well  documented  that  cycles  of  charge 
and  discharge  lead  to  capacity  fade  in  a  battery  over  its  lifetime 
[1-3].  Such  capacity  fade  is  suggested  to  be  partially  attributed 
to  mechanical  cracking  of  the  active  material  particles  induced  by 
repeated  lithium  intercalation/deintercalation  [4].  While  the  frac¬ 
ture  of  particles  has  been  observed  and  documented,  investigation 
of  processes  leading  to  crack  formation  and  development  of  models 
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capable  of  predicting  initiation  of  cracking  is  still  under  way  [5].  It 
should  be  mentioned  that  historically  the  area  of  battery  research 
was  driven  predominantly  by  electrochemical  expertise  and  mate¬ 
rials  development.  This  resulted  in  an  enormous  gap  in  research  on 
mechanical  behavior  of  electrode  materials.  At  the  present  stage 
it  is  evident  that  the  mechanical  engineering  side  of  the  problem 
needs  to  be  explored.  Determination  of  the  stress  state  of  a  particle 
during  charge/discharge  will  lead  to  understanding  the  origin  of 
fracture  and  will  provide  a  solution  in  order  to  improve  the  battery 
service  life. 

It  is  well  known  that  electrochemical  cycling  leads  to  capacity 
fade  in  lithium  ion  batteries,  which  creates  significant  limitations 
as  far  as  lifetime  and  reliability  of  a  product  is  concerned  [1,2].  In 
addition  to  processes  associated  with  chemical  reactions,  a  major 
source  of  capacity  fade  is  the  failure  of  electrode  particles  due  to  the 
stresses  developed  during  lithium  intercalation/deintercalation. 
Higher  lithium  insertion  rates  generally  result  in  higher  rate  of 
fracture  due  to  significant  lattice  strain  [5].  There  is  a  multitude 
of  experimental  evidence  supporting  the  notion  of  active  electrode 
particle  failure  due  to  intercalation  stresses  [6-9].  The  process  of 
cyclic  expansion  and  contraction  of  active  electrode  material  can 
be  viewed  from  the  standpoint  of  fatigue  and  damage  mechanics. 
However,  to  the  present  day,  the  term  “mechanical  fatigue”  has 
rarely  been  applied  to  the  processes  involved  in  electrochemical 
cycling. 

Involvement  of  mechanical  stresses  in  battery  degradation 
has  been  realized  only  recently  and  the  modeling  efforts  have 
started  utilizing  a  combination  of  diffusion  laws  with  continuum 
mechanics  [10-15].  Most  known  is  the  thermal  analogy  approach 
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Nomenclature 

c  lithium  concentration 

D  damage  parameter 

DLi  diffusion  coefficient 

e  dilatational  strain 

E  modulus  of  elasticity 

E+  modulus  of  elasticity  in  tension 

E~  modulus  of  elasticity  in  compression 

F  Faraday  constant 

h  material  constant  related  to  partial  closure  of  micro 

defects 

is  surface  current  density 

J  flux 

m  Weibull  modulus 

Pf!  probability  of  fracture 

r  radial  coordinate 

R  gas  constant 

R0  radius  of  the  particle 

t  time  variable 

T  temperature 

a  coefficient  of  thermal  expansion 

8y  Kronecker  delta 

£ij  components  of  strain  tensor 

er  radial  component  of  strain 

£t  tangential  component  of  strain 

r)k  k- th  eigenvector  of  stress  tensor 

v  Poisson  ratio 

<jSq  equivalent  stress 

oH  hydrostatic  stress 

cry  components  of  stress  tensor 

ok  k-th  eigenvalue  of  stress  tensor 

or  radial  component  of  stress 

crt  tangential  component  of  stress 

oVM  Von  Mises  stress 

crw  Weibull  stress 

{o)i  positive  part  of  stress  tensor 

{o)7.  negative  part  of  stress  tensor 

X  triaxiality  parameter 

£2  partial  molar  volume 

tr(or)  trace  of  tensor  a 

(•)  Macauley  bracket 


developed  by  Yang  [10]  for  calculation  of  chemical  stresses,  which 
has  been  adopted  in  subsequent  investigations  [11-13].  The  idea 
of  description  of  volumetric  expansion  due  to  diffusion  based  on 
equations  of  thermal  expansion  is  rather  appealing  due  to  simplic¬ 
ity  of  elastic  equations  and  this  approach  is  used  in  the  current 
work.  Computations  based  on  such  approach  can  give  a  good 
approximation  of  stress  distribution,  provided  an  appropriate  dif¬ 
fusion  law  is  available. 

The  electrochemical  cycling  of  the  electrochemical  cell  imposes 
cycles  of  expansion  and  contraction  similar  to  repeated  cycles  of 
thermal  straining.  Following  the  analogy,  the  notion  of  thermal 
fatigue  could  be  applied  to  failure  of  electrode  material.  The  phe¬ 
nomenon  of  thermal  fatigue  has  been  investigated  over  a  long 
period  of  time  [16-18]  and  the  presence  of  temperature  induced 
plastic  deformation  in  ductile  materials  has  been  reported  [16]. 
In  such  cases  the  relation  between  accumulation  of  damage  and 
plastic  strain  energy  results  in  the  fatigue  criterion  for  life  pre¬ 
diction  [19-21].  The  majority  of  Li-ion  battery  electrode  materials 
(Si,  LiMn204,  LiCo02,  representing  the  typical  anode  and  cathode 
compounds)  are  however  brittle  in  nature  and  the  accumulation  of 


damage  in  this  case  is  best  described  by  the  change  in  elastic  prop¬ 
erties  due  to  formation  of  microcracks.  When  the  microcracks  reach 
the  critical  size  (typically  very  small  compared  to  critical  cracks  in 
ductile  materials)  complete  failure  occurs.  In  this  regard  the  dam¬ 
age  can  be  related  to  the  fraction  of  Young’s  modulus  available  for 
further  elastic  deformation  and  the  value  of  modulus  equal  to  zero 
would  represent  complete  fracture.  Development  of  the  relation¬ 
ship  between  the  stress  state  and  damage  parameter  cannot  be 
based  on  Von  Mises  stress,  as  in  case  of  ductile  materials,  and  should 
include  a  different  formulation  of  equivalent  stress  [22,23]. 

The  present  work  represents  investigation  of  fracture  in  Si  par¬ 
ticles  used  as  an  active  material  for  anodes  in  a  Li-ion  battery. 
Without  venturing  into  details  of  electrochemical  aspects  of  bat¬ 
tery  operation,  the  problem  is  treated  from  the  stress  analysis  point 
of  view.  The  results  are  compared  with  the  acoustic  emission  data 
from  Si  electrode  cycling  and  a  good  correlation  was  found  between 
the  peak  of  emissions  and  the  predicted  time  when  the  damage 
parameter  reaches  its  maximum. 

2.  Evaluation  of  intercalation  induced  stress 

The  overall  problem  of  ionic  diffusion  in  solid  particles  is  very 
complex  from  both  electrochemical  and  mechanical  points  of  view 
and  therefore  necessary  approximations  are  made  to  reflect  partic¬ 
ular  aspects  of  the  overall  process.  By  considering  binary  diffusion 
coefficients  and  thermodynamic  driving  force  for  diffusion,  New¬ 
man  and  Christensen  [14,15]  developed  a  continuum  model  for 
computation  of  internal  stresses  in  a  spherical  particle  of  cathode 
material.  Another  approach  is  based  on  a  rather  natural  analogy 
between  volumetric  thermal  expansion  and  three-dimensional  dif¬ 
fusion  [10-13].  The  following  section  provides  brief  description  of 
the  thermal  analogy  model  in  terms  of  mathematical  formulation 
within  continuum  mechanics  framework.  The  particle  is  approxi¬ 
mated  as  a  homogenous  isotropic  sphere  for  simplicity.  The  process 
of  Li-ion  solid  state  diffusion  is  considered  as  a  process  analogous 
to  the  thermal  expansion  of  a  homogeneous  sphere  submerged 
into  liquid  media  which  has  higher  temperature  than  the  sphere. 
The  elastic  thermal  analysis  of  an  isotropic  media  is  governed  by 
Hooke’s  law  with  additional  thermal  expansion  term  [24] 

Eij  =  1  [(1  +  v)<jjj -vakk8jj]  +  aAT8i „  (1) 

Dummy  index  summation  convention  is  used  in  Eq.  (1 )  with  8y 
being  Kronecker  delta,  £y  and  Oy  represent  strain  and  stress  tensor 
components,  respectively,  v  is  the  Poisson  ratio,  a  is  the  coefficient 
of  thermal  expansion,  and  AT  represents  temperature  change  at  the 
material  point.  Analogous  to  Eq.  (1)  the  bulk  diffusion  in  isotropic 
media  can  be  represented  by  the  following  stress-strain  relation¬ 
ship  [10] 

Sij  =  | [(!  +  v)<7ij  -  vakkSij]  +  (2) 

where  £2  is  the  partial  molar  volume  and  c  is  the  concentration 
of  diffusing  species  at  the  material  point.  Since  £2  is  a  volumetric 
quantity  it  can  be  determined  analogously  to  volumetric  coefficient 
of  thermal  expansion  as 


e=T2Ac  (3) 

where  e  =  £kk  is  the  volumetric  strain. 

The  case  of  the  spherical  geometry  of  a  particle  produces  only 
two  components  of  the  stress  tensor  in  spherical  coordinates:  radial 
component  (crr)  and  tangential  component  ( ot ).  The  equilibrium 
condition  for  the  stress  tensor  in  case  of  a  spherical  particle  reduces 
to 

do r  2 .  /  a  \ 

+  -(^r  -Ot)  =  o  (4) 
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The  stresses  are  functions  of  elastic  strains  following  the  Hooke’s 
law 


or  = 


ot  = 


E 

(1  +  v)(l  -  2v) 
E 

(l  +  v)(l  -2v) 


(1  -  v)er  +  2vet-(l  +v) 
ver  +  £t-(l  +  v)  — 


£2c 

IT 


(5) 


where  the  components  of  strain  are  determined  from  the  radial 
displacement  as 


£r  — 


dU 

dr 


(6) 


Substitution  of  Eq.  (6)  in  expressions  for  stress  (5)  and  then  to 
equilibrium  Eq.  (4)  leads  to  the  ordinary  differential  equation  for 
displacement  as  a  function  of  radial  coordinate  r 


d2U  2  dU  2U  _  \  +  v  Q  dc 
dr2  +  r  dr  r2  ~  1  -  v  3  dr 


Provided  the  instantaneous  concentration  of  diffusing  species  is 
known,  solution  of  Eq.  (7)  with  back  substitution  in  Eqs.  (6)  and  (5) 
gives  the  stress  state  within  the  particle. 

The  diffusion  equation  is  introduced  as  an  extension  of  the  Fields 
law  with  the  component  corresponding  to  diffusion  driven  by  the 
hydrostatic  stress  oH  [10,11] 


=  DLiV2c  -  ^E(VcVaH  +  cV2aH ) 


(8) 


where  DLi  is  the  diffusion  coefficient,  R  is  the  gas  constant,  T  is  the 
absolute  temperature,  and  oH  =  okk/ 3  is  hydrostatic  stress.  When 
dependence  on  radial  coordinate  r  only  is  considered,  as  in  the  case 
of  isotropic  spherical  particle  with  uniform  distribution  of  solute  at 
the  surface,  the  diffusion  equation  becomes 


dc  _  I" d2c  2  9c  Q  9c  doH  Qc  f  d2oH  2  doH\ 
9t  =  u[dP+rW~Wdr^F~Wy~dP  +r~dFJ 


(9) 


Eqs.  (5)-(7)  together  with  the  diffusion  equation  (Eq.  (9)) 
describe  the  problem  of  ionic  diffusion  in  a  homogenous  spherical 
particle.  In  fact,  Eq.  (7),  being  an  ordinary  differential  equation,  can 
be  solved  analytically  which  results  in  the  following  expressions 
for  stress  components 


Or  = 


Ot  : 


2  QE 

3(1^5 

QE 


rR  o  rr 

1  /R30  /  cr2  dr  -  1  /r3  /  cr2  dr 
Jo  Jo 


3(1 -u) 


rKo  r 

2/Rq  /  cr2  dr  +  1  /r3  /  cr2  dr  -  c 

Jo  Jo 


(10) 


where  R0  is  the  radius  of  the  particle.  After  substituting  oH  =  ( or  + 
2ot)/3  with  stress  components  expressed  in  Eqs.  (9)  and  (10)  the 
differential  equation  for  concentration  evolution  becomes 


9c  (  92c  2  9c  /  9c 

at=Du  [dP  +  TTr+e[Tr 


n  ,'d2c  2  dc' 

+  9c{dP  +  TYr 


e  = 


2  n2E 
9RT(1  -  v) 


(ID 


=  const. 


Eq.  (1 1 )  represents  a  boundary-initial  value  problem  and  can  be 
solved  numerically  on  the  one-dimensional  space  grid.  It  should  be 
mentioned,  that  the  statement  of  the  problem  in  form  of  Eqs.  (10) 
and  (11)  is  performed  based  on  the  assumption  of  instantaneous 
mechanical  response,  i.e.  the  stress-strain  state  depends  on  time 
only  implicitly  via  concentration  function  c(r,  t).  At  each  particular 
time  step  during  electrochemical  cycling  process,  the  Eq.  (11)  is 
solved  first  and  the  stresses  are  obtained  as  a  function  of  radial 
coordinate  by  integration  of  Eq.  (10).  The  values  of  concentration 
of  Li  ions  from  the  previous  time  step  serve  as  initial  conditions 


Fig.  1.  Schematics  of  intercalation/deintercalation  stresses  in  a  spherical  particle. 


for  the  numerical  procedure  at  the  next  time  step.  The  boundary 
condition  imposed  at  the  surface  of  the  particle  is  controlled  by  the 
current  density  as  [1 1  ] 


J=  ~Du  (  Vc-^VaH 


r=R  o 


is 

E 


(12) 


where  F  is  the  Faraday  constant,  J  is  the  species  flux  and  is  is  the 
current  density  at  the  surface  of  the  particle.  Considering  solution 
for  stresses  (Eq.  (10))  this  boundary  condition  becomes 


y=  -dLi(i+0C)| 


r=R0 


}s 

F 


(13) 


The  second  boundary  condition  reflects  the  symmetry  of  con¬ 
centration  distribution  with  respect  to  the  center  of  the  particle 
and  is  represented  by  the  equation 


9c 

9r 


=  0 


r=0 


(14) 


3.  Damage  and  brittle  fracture 

Eq.  (11)  with  the  boundary  conditions  stated  in  Eq.  (13)  and 
Eq.  (14)  describes  the  evolution  of  concentration  of  Li  ions  in  a 
solid  homogenous  spherical  particle  with  time,  following  the  cyclic 
charge/discharge  curve.  Such  a  curve  is  typically  given  as  a  current 
(voltage)  versus  time  diagram  of  an  electrochemical  cell.  By  solv¬ 
ing  Eq.  (11)  the  concentration  as  a  function  of  a  radial  coordinate 
is  found  at  a  particular  time  and  the  stresses  are  obtained  by  inte¬ 
gration  of  expressions  in  Eq.  (10).  Thus  the  stress  distribution  is 
obtained. 

The  above  procedure  is  based  on  the  analogy  of  chemical  dif¬ 
fusion  stresses  to  the  thermal  stresses  in  continuum  mechanics.  In 
the  case  of  ductile  materials,  the  thermal  cyclic  loading  is  known 
to  contribute  to  the  fatigue  damage  accumulation  governed  by 
cyclic  plasticity  [16-18].  Many  of  the  typical  active  electrode  mate¬ 
rials  in  lithium  chemistries  are  brittle  in  nature,  such  as  LiCo02, 
LiMn204,  and  Si.  Therefore  the  damage  parameter  cannot  be  based 
on  accumulated  plastic  strain  energy  as  in  the  case  of  ductile 
engineering  materials  [21-23].  Electrode  material  deforms  volu- 
metrically  with  cyclic  insertion  and  removal  of  Li  ions  from  its 
crystallographic  structure.  Such  cycles  of  lithiation/delithiation 
produce  corresponding  compressive  and  tensile  stresses  within  a 
particle,  as  schematically  illustrated  in  Fig.  1. 

The  representation  of  strength  of  brittle  materials  (e.g.  ceramic) 
is  not  deterministic  in  nature  and  follows  the  Weibull  distribution 
[25,26].  The  assessment  of  tensile  strength  is  based  on  the  results 
of  load-to-rupture  tests  performed  on  large  set  of  specimens  with 
application  of  corresponding  statistical  analysis  [27].  The  likelihood 
of  fracture  is  governed  by  the  distribution  of  flaws  in  the  specimen 
and  is  based  on  the  weakest  link  approach.  The  distribution  of  flaws 
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is  assumed  to  be  stable  and  no  slow  crack  growth  is  occurring  dur¬ 
ing  the  loading  [27].  In  this  case  the  probability  of  fracture  of  a 
specimen  subject  to  tensile  stress  o  is  represented  as  [26] 

Pn  =  i  _  e-(<V<nv)m  (15) 

where  crw  and  m  are  the  Weibull  stress  and  the  Weibull  modu¬ 
lus  correspondingly  determined  from  the  series  of  load-to-rupture 
experiments. 

The  process  of  material  deformation  as  a  result  of  ionic  diffusion 
is  a  time-dependent  process  which  changes  material  itself  and  thus 
the  assumption  that  the  initial  population  of  internal  flaws  does  not 
change  with  time  is  no  longer  valid.  Formulation  for  damage  accu¬ 
mulation  as  the  material  is  subject  to  rising  stresses  due  to  species 
concentration  gradients  is  done  adopting  a  quasi-brittle  approach 
[26,28,29].  In  case  of  uniaxial  loading,  the  damage  D  is  described 
as  the  ratio  of  the  failed  area  of  the  sample  to  the  original  area.  If 
the  specimen  is  represented  as  an  assembly  of  N  identical  fibers 
out  of  which  n  have  failed  at  a  given  time  increment,  then  the  dam¬ 
age  is  n/N  and  represented  by  Eq.  (15).  The  effective  tensile  stress 
resulting  from  the  re-distribution  of  load  due  to  reduction  in  load- 
bearing  area  is  then  d  =  cr/(l  -  D).  The  effective  compressive  stress 
is  cr  =  cr/(l  -  hD)  and  the  constant  h  reflects  the  partial  recovery  of 
elasticity  in  compression  due  to  closure  of  micro-cracks.  The  effec¬ 
tive  Young’s  modulus  in  tension  and  compression  is  represented 
as 

/  E+=E(1-D) 

\E-=E(l-hD)  (16) 

and  the  constant  h  is  therefore  h  =  (E  -E~)/{E  -£+).  For  most 
materials,  the  ratio  h  has  a  value  close  to  0.2  [26]  and  this  value 
is  adopted  in  the  present  investigation. 

The  above  reasoning  can  be  extended  to  the  case  of  multiaxial 
stress  state  in  which  case  the  effective  stress  tensor  is  introduced. 
The  positive  and  negative  parts  of  a  stress  tensor  are  constructed 
from  its  positive  (negative)  eigenvalues  ak  and  corresponding 
eigenvectors  fjk  as  follows 


Fig.  2.  Si  particles  used  in  the  experiments. 


where  ovm  is  the  Von  Mises  stress  and  x  is  the  triaxiality  parameter. 
With  inclusion  of  elasticity  change  due  to  damage  the  equivalent 
stress  is  expresses  by 

creq  =  yj{y  +v)(CT)J(cr}J-  v(trcr>2  +  [(1  +  v)<ct>-(ct>-  -  v(-trcr}2]  (21) 

It  can  be  seen  that  the  damage  evolution  and  stress  state  are 
interrelated  by  Eqs.  (19)  and  (21 )  via  the  equivalent  stress  parame¬ 
ter.  The  overall  solution  procedure  can  be  described  as  follows.  The 
electrochemical  “loading  curve”  is  represented  by  the  variation  of 
current  with  time.  At  the  initial  moment  the  initial  concentration  of 
Li  inside  the  particle  is  known  and  the  stresses  can  be  obtained  from 
Eq.  (10).  The  initial  damage  is  taken  to  be  zero  and  then  corrected 
by  Eq.  (19)  using  the  equivalent  stress  from  Eq.  (21 ).  The  solution  is 
continued  following  the  same  steps  for  each  point  along  the  current 
versus  time  curve.  Following  the  assumption  of  isotropic  mate¬ 
rial,  the  damage  is  isotropic  as  well  and  represented  by  one  scalar 
variable. 


H  (a*) 


k= 1 


(a)7  =  Vy  -  <u)+ 


(17) 


where  ()  represents  Macauley  brackets  such  that  (x)  = 
f  x,  ifx  >  0 
\0,  ifx  =  0‘ 

The  above  definition  results  in  the  following  expression  for  the 
effective  stress  tensor  in  multiaxial  stress  state 


-  + 


(O') 


1  -D  1  -hD 


-H~( 

1  -2v  V 


8ki(v)u  ~  (<*kk)  8kiWkl  +  (-cruc)  \ 

1  -D  +  1  -  hD  J 


Sij  (18) 


As  for  the  damage  parameter  under  multiaxial  stress,  the  quan¬ 
tity  a  in  Eq.  (15)  should  be  substituted  with  a  corresponding 
equivalent  stress  parameter  representing  the  current  stress  state. 
Thus  the  damage  parameter  becomes  represented  by 


D=1  -  e-(°eq/ow)m 


(19) 


We  choose  the  equivalent  stress  definition  in  form  that  depends 
on  the  hydrostatic  stress  [26] 


&eq  =  &VMVX 

X  =  |(l+v)  +  3(l-2 

3  \VvmJ 


4.  Comparison  with  experimental  data 

In  the  present  section  the  computational  results  of  the  dam¬ 
age  accumulation  are  compared  to  the  results  of  acoustic  emission 
(AE)  experiments  on  anodes  with  Si  active  material.  A  series  of 
experiments  were  conducted  on  lithium  half  cells  with  Li  counter 
electrode  under  different  charge/discharge  conditions.  The  sili¬ 
con  particles  (-100/+325  mesh  from  Alfa  Aesar)  were  mixed  with 
10wt%  PVDF  binder  and  10wt%  carbon  black  as  conductive  addi¬ 
tive,  and  spread  over  copper  circular  disks  with  radius  equal  to 
12  mm.  The  surface  density  of  Si  was  2.5  mg/cm2.  Silicon  particles 
in  as  received  condition  are  shown  on  the  SEM  picture  in  Fig.  2.  The 
average  size  of  the  particles  was  100  p,m. 

The  half  cell  employed  in  the  experiments  consisted  of  two 
stainless  steel  plates  and  insulating  PTFE  gasket  in  which  the  Si 
composite  anode,  separator  and  Li  foil  were  placed.  AIM  solution 
of  LiPF6  in  ethylene  carbonate  and  ethyl  methyl  carbonate  (3:7  wt) 
(Ferro  Corp.)  was  used  as  electrolyte.  The  half  cells  were  cycled 
using  Maccor  battery  tester.  The  AE  data  was  collected  using  a  4 
channel  system  from  Physical  Acoustics. 

The  experiments  were  conducted  under  constant 
current-constant  voltage  (CCCV)  conditions  with  correspond¬ 
ing  cut-off  values.  CCCV  process  consists  of  periods  of  time  when 
the  cell  is  held  at  constant  current  until  the  voltage  reaches  desired 
cut-off  value  and  then  the  cell  is  held  at  constant  voltage  until  the 
current  drops  to  the  designed  value.  The  details  of  CCCV  procedure 
can  be  found  in  [30,31  ].  The  representative  cyclic  curve  is  shown  in 
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Fig.  3.  Example  of  constant  current-constant  voltage  cycling  curve  used  in  the 
experiments. 


Fig.  3.  Total  of  eleven  cells  with  Si  anodes  were  cycled  with  current 
amplitudes  ranging  from  0.29  to  2.3  mA. 

In  order  to  eliminate  noise,  acoustic  emissions  representing 
particle  cracking  were  filtered  based  on  21  dB  threshold.  The 
electromagnetic  interference  (EMI)  events  were  registered  by  an 
additional  AE  sensor  placed  near  the  cell  and  these  events  were 
filtered  out.  The  acoustic  emission  parameters  (i.e.  duration,  ampli¬ 
tude  and  frequency)  were  compared  to  those  available  in  the 
literature  [32-34]  and  a  good  agreement  was  observed.  The  details 
of  signal  collection  and  post-processing  can  be  found  in  corre¬ 
sponding  Refs.  [31-34]. 

Determination  of  several  material  constants  is  required  in  order 
to  use  the  equations  specified  above.  The  partial  molar  volume 
was  determined  based  on  the  analogy  with  volumetric  coefficient 
of  thermal  expansion.  The  volumetric  strain  (ej)  due  to  the  tem¬ 
perature  change  AT  is  expressed  as  eT  =  av  AT,  where  av  is  the 
volumetric  coefficient  of  thermal  expansion.  Analogously,  for  volu¬ 
metric  expansion  due  to  intercalation  of  lithium,  ec  =  £2  AC,  where 
AC  is  the  change  in  lithium  concentration.  Based  on  300%  volume 
increase  in  Li*Si  when  x  changes  from  0  to  4.4  [30]  the  value  of  Q 
is  7.69  x  10-6  m3  mol-1.  In  the  current  work  the  diffusion  coeffi¬ 
cient  is  assumed  to  be  constant  and  independent  of  state  of  charge 
(stress  state).  The  parameters  of  Weibull  distribution  for  strength 
of  Si  were  adopted  from  the  literature  [25].  The  apparent  value 
of  the  diffusion  coefficient  and  other  material  constants  is  shown 
in  Table  1  together  with  corresponding  references  reporting  these 
values. 

The  solution  of  the  diffusion  equation  (Eq.  (11))  is  controlled 
by  the  surface  current  density  boundary  condition  (Eq.  (13)).  At 
each  point  in  time,  the  current  density  at  the  particle  surface  was 
obtained  by  dividing  the  total  current  by  the  average  number  of 
Si  particles  in  the  electrode  and  by  the  average  surface  area  of  a 
particle.  Thus  each  point  from  the  current  cycling  curve  (Fig.  3)  can 
be  converted  to  an  average  surface  current  density  of  a  particle. 
Each  time  interval  from  the  cycling  electrical  data  was  split  into 
sub-intervals  for  time  integration  of  Eq.  (11).  The  numerical  solu¬ 
tion  of  initial-boundary  value  problem  was  performed  by  backward 
time  centered  space  (BTCS)  method.  The  overall  arrangement  of  the 
solution  steps  is  shown  in  Table  2,  where  indexes  k  and  l  indicate 
time  step  and  node  numbers  respectively.  The  primary  interests 
were  the  evolution  of  maximum  damage  within  a  particle  with 
time,  distribution  of  stress  and  distribution  of  damage  inside  the 


Table  1 

Material  properties  used  in  calculations. 


Du  m2  s-1 

Q,  m3  mol-1 

EG  Pa 

V 

crw  GPa 

m 

p  kg  rrr 3 

5.1  x  lCH6 

7.69  x  10“6 

190 

0.218 

6.1 

3.6 

2.33  x10s 

[36-39] 

[30] 

[25,35] 

[35] 

[25] 

[25] 

[35] 

Table  2 

Solution  steps  for  determination  of  damage  parameter. 

1  Input  cyclic  curve  Ik,  tk 

2  t  =  0,  Do  =  0,  Cqi  =  0 

3  For  each  tk : 

3.1  Obtain  4k)  for  BC  Eq.(13) 

3.2  Solve  FDE  form  of  Eq.  (11)  for  nodal  cki  (Thomas  algorithm  is  used  for 
tridiagonal  matrix). 

3.3  £(*’') -Eq.  (10) 

3.4<r£'>-Eq.(21) 

3.5DU-Eq.(19) 

3.6  Search  for  maxDk 

3.7  t  =  tk  +  At  and  go  to  3.1 


Fig.  4.  AE  data  for  electrochemical  cycling  of  Si  anode:  specimen  ID  48. 


particle.  Influence  of  cycling  conditions  on  damage  development 
was  investigated.  The  results  are  compared  with  experimental  AE 
data  records  of  cracking  events. 

As  an  example,  the  acoustic  emission  data  for  a  representative 
specimen  (sample  ID  48)  is  shown  in  Fig.  4  together  with  the  cyclic 
current  as  a  function  of  time.  The  markers  represent  the  ampli¬ 
tude  of  the  acoustic  events  (hits)  distributed  within  the  duration 
of  experiment.  It  can  be  seen  that  most  of  the  cracking  associated 
with  the  AE  hits  occurs  during  the  first  discharge,  i.e.  at  the  begin¬ 
ning  of  the  first  cycle.  This  situation  appears  to  be  typical  for  all  of 
the  specimens  investigated.  The  subsequent  emissions  registered 
at  the  beginning  of  the  corresponding  charge  (positive  current)  and 
discharge  (negative  current)  portions  of  the  cyclic  curve  can  be 
attributed  to  the  formation  of  new  cracks  within  the  grains  frac¬ 
tured  during  the  first  cycle.  The  influence  of  the  newly  created 
cracks  on  the  electrochemical  and  mechanical  properties  of  the 
particle  is  beyond  the  scope  of  current  investigation.  At  the  present 
stage  it  appears  to  be  essential  to  verify  that  the  approach  described 
above  can  predict  the  fracture  of  particle  in  agreement  with  the 
acoustic  emission  data  during  the  first  cycle.  That  is  the  “safe-life” 
approach  is  taken  and  only  the  accumulation  of  damage  in  the  ini¬ 
tially  not  fractured  particle  is  considered.  Under  this  condition,  the 
first  maximum  of  emission  rate  should  coincide  with  the  maximum 
of  the  damage  parameter. 

The  damage  evolution  with  time  is  shown  together  with  the  rate 
of  emissions  in  Fig.  5  for  the  specimen  ID  48.  Only  the  first  discharge 
portion  of  the  cycle  is  shown.  Since  the  silicon  electrode  is  cycled 
against  lithium,  the  former  becomes  a  cathode  in  the  half-cell  and 
the  intercalation  of  Li  ions  occurs  during  discharge.  The  dash-dotted 
line  in  Fig.  5  represents  the  rate  of  the  emissions  collected  during 
the  experiment  and  the  solid  line  represents  the  calculated  damage 
based  on  the  current  evolution  in  time.  Emission  rate  was  calculated 
by  differentiating  the  number  of  events  versus  time  data.  The  data 
is  most  successfully  fitted  with  the  Hill  equation 


F{x)  =  b  + 


m-b 

1  +(V*)r 


(22) 
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Fig.  5.  Correlation  between  emissions  rate  from  experiment  and  computed  evolu¬ 
tion  of  maximum  damage. 

where  b ,  m,  Ih,  and  r  are  fitting  parameters.  The  point  where  the 
rate  of  acoustic  emission  events  reaches  its  maximum  is  consid¬ 
ered  to  be  the  point  where  most  of  the  particles  within  electrode 
experience  fracture. 

With  the  use  of  the  damage  parameter,  the  time  when  D  =  1  rep¬ 
resents  the  time  to  failure  of  a  particle.  It  can  be  seen  from  Fig.  5 
that  the  maximum  of  damage  parameter  (value  of  1.0)  correlates 
well  with  the  maximum  of  emissions  rate. 

The  correlation  between  experimental  data  and  predicted  frac¬ 
ture  based  on  D  =  1  criterion  is  shown  in  Fig.  6.  The  markers 
represent  data  for  individual  specimens.  The  solid  line  indicates 
perfect  correlation  and  dashed  lines  show  the  boundary  within  the 
factor-of-two  dispersion  from  the  perfect  correlation.  It  can  be  seen 
that  the  experimental  data  is  in  good  agreement  with  computed 
values  of  time  to  failure  based  on  the  brittle  damage  parameter. 

Out  of  eleven  specimens  under  investigation,  ten  were  cycled 
between  1.3  V  and  50  mV  under  different  cut-off  currents.  One 
specimen  was  tested  with  higher  cut-off  voltage,  and  cycled 
between  1.3  and  170  mV.  The  particles  in  this  electrode  did  not 
fracture  and  no  acoustic  emission  was  detected  during  cycling. 
This  specimen  is  indicated  by  an  arrow  in  Fig.  6  implying  that  the 
electrode  particles  did  not  fail.  An  arbitrary  time  of  80  x  104  s  was 


Fig.  6.  Correlation  between  observed  and  predicted  fracture  times;  an  arrow  next 
to  the  data  point  indicates  that  no  fracture  was  observed. 


Fig.  7.  Computed  damage  evolution  as  a  function  of  time  and  radial  coordinate  of  a 
spherical  particle. 


chosen  to  represent  the  infinite  time  to  fracture.  Fligher  cut-off 
voltage  resulted  in  only  partial  lithiation  of  Si  which  reduced  the 
internal  stresses  within  particles  thus  reducing  the  damage.  The 
maximum  damage  in  this  specimen  reached  the  value  of  0.15  only, 
therefore  predicting  the  experimentally  observed  non-failure. 

In  Fig.  5  the  maximum  value  of  damage  within  the  particle  is 
shown  for  each  corresponding  moment  of  time.  The  distribution  of 
damage  within  the  particle  is  of  particular  interest  since  the  loca¬ 
tion  of  its  maximum  would  represent  the  point  where  the  fatal 
crack  originates.  The  damage  as  a  function  of  the  radial  coordi¬ 
nate  and  time  (D(r,t))  is  plotted  in  Fig.  7.  The  three-dimensional 
representation  in  Fig.  7  depicts  the  evolution  of  damage  profiles 
with  time  as  the  lithiation  proceeds  following  the  initial  discharge 
curve  in  Fig.  3.  The  particle  size  was  chosen  as  40  |xm  and  thus 
the  radial  coordinate  changes  from  zero  (center  of  the  particle) 
to  20  fjim  (surface  of  the  particle).  Since  the  damage  at  the  initial 
stages  of  discharge  lithiation  is  close  to  zero,  the  portion  of  the 
damage  surface  corresponding  to  onset  of  lithiation  is  not  plotted. 
The  intersection  of  the  surface  with  the  damage-radial  coordinate 
planes  represents  the  two  dimensional  plots  of  the  function  D(r) 
corresponding  to  different  points  in  time.  The  intersection  of  the 
surface  with  damage-time  plane  where  damage  attains  its  maxi¬ 
mum  represents  time  evolution  Dmax(t)  similar  to  that  shown  in 
Fig.  5.  The  computation  continued  until  the  value  of  D(r,t)  reached 
its  maximum  and  thus  the  damage  profile  at  the  last  point  in  time 
represents  the  spatial  distribution  of  damage  at  the  moment  of  frac¬ 
ture.  It  can  be  seen  that  within  the  entire  time  span,  the  maximum 
damage  occurs  at  the  center  of  the  particle  which  is  where  the  peak 
of  damage  is  situated  when  the  cracking  initiates.  Thus  the  particle 
is  expected  to  fracture  in  the  center  when  the  lithium  is  inserted. 

The  results  reported  above  are  based  on  the  experimental  data 
suggesting  that  the  peak  of  acoustic  emissions  occurs  during  the 
first  discharge;  that  is  most  of  the  particles  fracture  during  the 
initial  insertion  of  Li.  The  calculations  show  rapid  increase  of  the 
damage  and  the  points  of  D  =  1  coincide  well  with  the  peaks  of 
acoustic  emission  (Fig.  6).  The  predictions  are  based  on  the  aver¬ 
age  radius  of  the  particle  (50  [xm)  used  in  the  experiments.  At  this 
point  it  seems  only  logical  to  investigate  the  effect  of  particle  size 
on  damage  evolution  with  time.  Once  again  the  loading  condi¬ 
tions  applied  to  specimen  48  were  used  as  the  input  values.  The 
damage-time  curves  for  the  particles  of  three  different  sizes  are 
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Fig.  8.  Distribution  of  maximum  damage  for  three  different  radii  of  Si  particle. 

shown  in  Fig.  8.  The  time  period  shown  in  Fig.  8  corresponds  to  the 
initial  discharge  part  of  the  loading  curve  (Fig.  3).  Two  trends  can  be 
clearly  seen  from  Fig.  8.  The  value  of  maximum  damage  decreases 
with  decrease  in  particle  radius  and  the  time  to  the  peak  of  the 
damage  also  decreases  as  the  particle  size  is  reduced.  The  former 
observation  is  rather  intuitive  as  the  larger  particle  would  produce 
greater  stress  fields  which  would  lead  to  fracture.  The  location  of 
the  damage  maxima  closer  to  the  onset  of  discharge  follows  from 
the  fact  that  the  shorter  time  is  needed  for  lithiation  of  a  smaller 
particle. 

The  continued  investigation  of  the  size  effect  is  shown  in  Fig.  9, 
where  the  maximum  damage  and  the  time  until  this  maximum 
damage  is  reached  is  plotted  against  the  particle  radius  ranging 
from  1  to  50  p,m.  The  vertical  dashed  line  represents  boundary 
between  two  particle  sizes  identified  by  the  mesh  numbers  and 
will  be  explained  in  the  Discussion  section  of  the  paper.  It  can  be 
seen  that  as  the  particle  size  decreases  so  does  the  damage,  reach¬ 
ing  the  limit  close  to  zero  for  infinitely  small  particle.  The  damage 
then  increases  rapidly  within  the  range  of  particle  radii  between  1 0 
and  20  p,m  and  asymptotically  reaches  the  plateau  with  the  value 
of  1.0.  Finally,  it  can  be  seen  that  for  particles  smaller  than  10  [xm 
in  radius  the  value  of  damage  is  rather  negligible  suggesting  that 
the  particles  of  this  range  of  sizes  will  not  fracture  upon  the  ini¬ 
tial  discharge.  It  should  be  mentioned  that  the  scale  of  10  pun  is  far 
greater  than  nanoscale  typically  considered  to  be  “stress  safe”  on  a 
single  particle  level  [2,3]. 

5.  Discussion 

The  present  investigation  represents  the  analysis  of  diffusion 
induced  fracture  of  Si  particles  during  lithium  insertion  and  extrac¬ 
tion  based  on  a  working  principle  of  a  Li-ion  battery.  The  problem 
is  treated  from  the  continuum  mechanics  point  of  view.  The  prob- 


Fig.  9.  Damage  parameter  as  a  function  of  particle  size. 


Fig.  10.  SEM  of  a  cycled  Si  particle. 


abilistic  formalism  of  deformation  analysis  of  brittle  materials 
resulted  in  damage  parameter  reflecting  the  probability  of  frac¬ 
ture  with  D  ->  1  representing  the  prediction  of  complete  failure. 
Within  this  notion  the  stress  state  is  in  turn  affected  by  the  change 
in  the  elastic  modulus  due  to  damage  following  Eq.  (19).  Distinction 
between  tension  and  compression  is  made  through  consideration 
of  microcracks  closure  when  the  stress  state  is  compressive. 

Examination  of  AE  data  from  electrodes  tested  (Fig.  4)  reveals 
the  maximum  density  of  emission  events  during  the  first  discharge 
of  the  cyclic  curve.  In  a  half-cell,  an  anode  formally  serves  as  a  cath¬ 
ode  as  it  is  cycled  against  lithium  metal  and  thus  discharge  portion 
of  a  cycle  represents  insertion  of  Li  ions  into  Si  particles.  While 
the  highest  density  of  emissions  occurs  during  the  first  discharge, 
Fig.  4  shows  spikes  of  emissions  during  subsequent  cycling.  The  AE 
events  occur  during  both  the  discharge,  when  the  stress  state  on 
the  surface  is  compressive  as  well  as  during  charging,  when  the 
particle  surface  stress  state  is  tensile.  At  the  present  moment  it  is 
rather  difficult  to  make  a  definite  statement  as  for  the  cause  of  sub¬ 
sequent  bursts  of  emissions.  It  could  be  suggested  that  the  initial 
burst  of  emissions  is  due  to  cracking  while  later  on  in  subsequent 
cycles  the  fretting  of  the  existing  fracture  surfaces  and  small  incre¬ 
ments  of  crack  growth  can  be  attributed.  Since  the  examination 
of  electrode  material  is  done  ex  situ  it  is  impossible  to  determine 
the  exact  sequence  of  events  leading  to  final  fracture  of  a  particle 
into  pieces.  The  continuation  of  the  current  research  will  be  done 
by  in  situ  optical  microscopy  using  a  specially  designed  Li  half  cell. 
Such  in  situ  observations  will  give  valuable  information  regarding 
development  of  cracking  with  electrochemical  cycling. 

An  example  of  cycled  particle  is  shown  on  SEM  image  in  Fig.  10. 
Multiple  cracks  separating  the  particle  into  smaller  blocks  can 
be  observed,  however  the  location  of  crack  initiation  cannot  be 
deducted  from  this  picture.  It  is  quite  possible  that  the  cracking 
during  the  initial  discharge  started  from  the  center  of  the  particle 
where  the  prediction  places  the  maximum  of  the  damage  parame¬ 
ter  (Fig.  7)  separating  the  particle  into  several  blocks.  Subsequent 
o  charge  would  introduce  the  surface  cracks  in  those  blocks  produc¬ 
ing  the  burst  of  emissions  corresponding  to  the  charging  part  of 
a  the  cycle.  Then  during  the  next  discharge  the  internal  cracks  would 
appear  again  and  so  on,  until  the  block  reaches  a  critically  small 
size  to  develop  enough  stress  for  fracture.  The  geometry  of  the 
internal  crack  produced  during  the  first  discharge  depends  on  the 
distribution  of  initial  internal  defects  in  the  particle  and  is  impossi¬ 
ble  to  predict  within  the  generalized  homogeneous  spherical  body 
approach.  In  other  words  it  is  impossible  to  quantify  the  amount  of 
blocks  produced  by  initial  fracture  of  a  particle. 
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Fig.  11.  Comparison  between  AE  cumulative  events  from  two  different  particle  size 
distributions. 
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Investigation  of  the  particle  size  effect  on  the  evolution  of  dam¬ 
age  produces  the  result  which  reveals  the  existence  of  particle 
size  range,  for  which  the  value  of  maximum  damage  is  less  than 
1  (Fig.  9).  Additional  examination  of  Fig.  9  yields  observation  that 
as  the  particle  radius  comes  close  than  10  pan  the  damage  becomes 
negligibly  small.  The  particle  size  in  the  vicinity  of  10  jxm  in  radius 
thus  can  be  considered  as  fracture  safe  under  the  current  cycling 
conditions. 

It  should  be  mentioned  that  the  above  size  range  is  still  much 
larger  than  expected  in  the  electrochemical  community  where  the 
shift  towards  nano-sized  electrode  particles  is  considered  to  be  one 
of  the  solutions  for  extension  of  battery  life  [2,3,40].  The  benefits 
of  nanoscopic  materials  for  lithium  battery  applications  consist  of 
short  diffusion  distances,  which  implies  increase  in  rate  capabil¬ 
ity,  and  large  surface  areas,  which  allows  for  fast  absorption  of 
lithium  ions  [3].  Less  severe  concentration  gradients  inside  smaller 
particles  result  in  reduced  magnitude  of  stress  which  prevents  frac¬ 
ture.  It  has  been  reported  that  reduction  of  LiMn1.5Nio.5O4  particle 
size  from  1  p,m  to  70  nm  resulted  in  significant  improvement  in 
capacity  retention  [41].  Work  on  nano-structured  arrays  of  Si  pil¬ 
lars  (~600nm  in  diameter)  showed  that  at  this  length  scale  the 
material  maintained  its  structural  integrity  during  electrochemical 
cycling  [42].  The  results  of  current  investigation  suggest  that  pre¬ 
vention  of  cracking  may  be  achieved  in  particles  with  size  much 
larger  than  the  nanometer  region. 

In  order  to  validate  the  size  effect  shown  in  Fig.  9,  one  of  the 
electrodes  was  prepared  with  Si  particles  of  -325  mesh,  that  is  less 
than  22  |xm  in  radius.  The  electrode  was  cycled  between  1.3  V  and 
50  mV  with  0.53  mA  current  amplitude.  The  results  are  shown  in 
Fig.  11  and  compared  to  the  AE  results  from  the  electrode  with 
-100/+325mesh  particles  (average  radius  50|jim)  cycled  under 
the  same  conditions.  According  to  ASTM-E1 1  standard,  the  mesh 
numbers  are  used  in  description  of  sizes  of  particulate  materials 
and  denote  the  sieve  openings  through  which  particles  are  passed. 
The  sign  in  front  of  the  number  corresponds  to  whether  particles 
go  through  the  sieve  or  are  trapped.  In  our  case,  the  designation 
-325  mesh  indicates  that  particles  smaller  than  44  pan  will  pass 
through  the  sieve,  and  designation  -100/ +325  means  that  parti¬ 
cles  will  pass  through  1 50  p,m  openings  ( - 1 00)  but  will  be  stopped 
by  44fxm  sieve  (+325).  The  line  corresponding  to  particle  diam¬ 
eter  of  44  |jim  and  thus  demarcating  the  two  meshes  is  placed  in 


Fig.  9  in  order  to  distinguish  the  corresponding  values  of  damage. 
It  can  be  seen  in  Fig.  9  that  the  part  of  the  maximum  damage  curve 
corresponding  to  the  particles  smaller  than  44  p,m  is  distinguished 
by  the  abrupt  drop  in  values  of  damage  parameter,  while  for  the 
particles  larger  than  44  p,m  the  damage  approaches  the  asymp¬ 
totic  maximum  value  of  unity.  Confirming  the  predictions  from 
Fig.  9,  Fig.  1 1  shows  that  the  number  of  hits  registered  in  the  small- 
particle  electrode  is  magnitudes  less  than  in  the  case  of  50  p,m 
particles.  The  experimental  work,  which  targets  investigation  of 
size  effect  on  fracture  time,  is  under  way  and  should  further  reveal 
the  dependence  of  damage  on  particle  size  as  well  as  demonstrate 
the  existence  of  critical  size  following  the  computational  results 
presented  in  current  manuscript. 


6.  Conclusions 

The  problem  of  Si  particles  fracture  upon  insertion  of  lithium 
in  a  Li-ion  half-cell  was  treated  within  the  continuum  mechan¬ 
ics  approach  with  brittle  damage  parameter.  The  results  predict 
fracture  upon  the  immediate  first  discharge  and  are  in  agree¬ 
ment  with  experimental  acoustic  emissions  data  obtained  from 
charge/discharge  cycling  of  a  half  cell  with  Si  as  anode  active  mate¬ 
rial.  Consideration  of  particle  size  effect  on  developed  chemical 
stresses  suggests  a  critical  size  of  a  particle  which  would  not  frac¬ 
ture  upon  initial  lithiation.  Such  critical  size  appears  to  be  within 
micrometer  scale. 
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